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Abstract 

A Wireless Energy Harvesting Node (WEHN) operating in linear vector Gaussian channels with arbitrarily 
distributed input symbols is considered in this paper. The precoding strategy that maximizes the mutual information 
along N independent channel accesses is studied under non-causal knowledge of the channel state and harvested 
energy (commonly known as offline approach). It is shown that, at each channel use, the left singular vectors of the 
precoder are equal to the eigenvectors of the Gram channel matrix. Additionally, an expression that relates the optimal 
singular values of the precoder with the energy harvesting profile through the Minimum Mean-Square Error (MMSE) 
matrix is obtained. Then, the specific situation in which the right singular vectors of the precoder are set to the identity 
matrix is considered. In this scenario, the optimal offline power allocation, named Mercury Water-Flowing, is derived 
and an intuitive graphical representation is presented. Two optimal offline algorithms to compute the Mercury Water- 
Flowing solution are proposed and an exhaustive study of their computational complexity is performed. Moreover, an 
online algorithm is designed, which only uses causal knowledge of the harvested energy and channel state. Finally, 
the achieved mutual information is evaluated through simulation. 
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I. Introduction 



Battery powered devices are becoming broadly used due to the high mobility and flexibility provided to users. 
As Moore predicted in 1965 [1], nodes' processing capability keeps increasing as transistors shrink year after year. 
However, the growth of battery capacity is slower and thus energy availability is becoming the bottleneck in the 
computational capabilities of wireless nodes. 

Energy harvesting, which is known as the process of collecting energy from the environment by different means 
(e.g. solar cells, piezoelectric generators, etc.), has become a potential technology to charge batteries and, therefore, 
expand the lifetime of battery powered devices (e.g., handheld devices or sensor nodes), which we refer to as 
WEHNs. The presence of energy harvesters implies a loss of optimality of the traditional transmission policies, 
such as the well-known Waterfilling (WF) strategy [2], because the common transmission power constraint must 
be replaced by a set of Energy Causality Constraints (ECCs), which impose that energy must be harvested before 
it can be used by the node. 

In general, the energy harvesting process is modeled as a set of energy packets arriving to the node at different 
time instants and with different amounts of energy. 1 There exist two well established approaches for the design of 
optimal transmission strategies, namely, online and offline. The online approach assumes that the node only has 
some statistical knowledge of the dynamics of the energy harvesting process, which can be realistic in practice. The 
offline approach assumes that the node has full knowledge of the amount and arrival time of each energy packet, 
which is an idealistic situation that provides analytical and intuitive solutions and, therefore, it is a good first step 
to gain insight for the later design of the online transmission strategy. Using this model for the energy harvesting 
process, references [3]— [10] derived the optimal resource allocation in different scenarios. The authors of [6] and [7] 
found the power allocation strategy, named Directional Water-Filling (DWF), that maximizes the total throughput 
of a WEHN operating in a point-to-point link, i.e., 



where n is the channel use index, A„ is the channel gain, Tj is a set that contains the channel accesses between 
two consecutive energy arrivals, and Wj is the water level associated to Tj. 2 

For nodes without energy harvesting capabilities, the capacity of linear vector Gaussian channels was given 
in [12], where it was shown that given a certain constraint in the transmitted power, Pc, capacity is achieved by 

'Observe that any energy harvesting profile can be accurately modeled by making the inter arrival times sufficiently small. 

2 Notation: Matrices and vectors are denoted by upper and lower case bold letters, respectively. [v] n denotes the n-th component of the vector 
v. [A]p<j is the component in the p-th row and q-th column of matrix A. I n is the identity matrix of order n, l n is a column vector of n 
ones. DxF denotes the Jacobian of the matrix function F with respect to (w.r.t.) the matrix variable X [11]. The superscript ( ) T denotes the 
transpose operator. diag(X) is the column vector that contains the diagonal elements of the matrix A. Diag(v) is a diagonal matrix where the 
entries of the diagonal are given by the vector v. vec(X) returns a vector that stacks the columns of X. Tr(-) denotes the trace of a matrix. 
var{ } is the variance of some random variable. The Kronecker product is denoted by the symbol <X). Finally, = max{0, i}. 
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diagonalizing the observed channel in K independent streams. Then, a Gaussian distributed codeword is transmitted 
over each stream whose power is obtained from the well known WF solution: 



The main difference between DWF (1) and WF (2) is that in the former the water level depends on the channel 
access under consideration. 

Both (1) and (2) are optimal when the distribution of the input is Gaussian. However, in practical scenarios, finite 
constellations are used instead of the ideal Gaussian signaling, e.g., Q-PAM and Q-QAM, where Q denotes the 
alphabet cardinality. In the low Signal to Noise Ratio (SNR) regime, the capacity achieved with finite constellations 
is very close to the one achieved by Gaussian signaling. However, the mutual information asymptotically saturates 
when the SNR increases as not more than log 2 Q bit per channel use can be sent (see Fig. 1 in [13]). This must 
be taken into account in the design of the optimal power allocation when the input symbols are constrained to 
belong to a finite alphabet. In opposition to the Gaussian case, where the better the channel gain, the higher the 
allocated power, when arbitrary constellations are used, there exists a tradeoff between the alphabet cardinality and 
the channel gain. In [14], the optimal power allocation was found for a node without energy harvesting capabilities 
and with arbitrary distributed input symbols. To do so, the authors of [14] used the relation between the mutual 
information and the MMSE, which was revealed in [15] and further generalized in [16], as summarized in the 
following lines. 

In [15], Guo et al. revealed that the derivative of the mutual information with respect to (w.r.t.) the SNR for a 
real-valued scalar Gaussian channel is proportional to the MMSE, i.e., 



where x is the channel input, n is the observed noise and mmse(snr) = E{(ir— x) 2 }, where x = K{x\ ^snrx+n} is 
the conditional mean estimator. The mutual information in linear vector Gaussian channels was further characterized 
in [16], where its partial derivatives w.r.t. arbitrary system parameters were determined, e.g., the gradient w.r.t. the 
channel matrix, H, was found to be 



where x is the vector input, n is the noise, E = E{(x — x)(x — x) T } is the MMSE matrix, and x = E{x|Hx + n}. 

Thanks to the relationship in (3), the power allocation that maximizes the mutual information over a set of parallel 
channels (each of them denoted by a different index k) with finite alphabet inputs was derived in [14] and named 




(2) 




(3) 
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Mercury/Waterfilling (H g WF), i.e., 

where Gfc(^) is the mercury factor that depends on the input distribution and is defined as 

fl-mmje^W if [0,1], 
G fc (^) = < (6) 
[l if -0 > 1- 

This result showed that the optimal power allocation not only depends on the channel gain as in the Gaussian 
signaling case, but also on the shape and size of the constellation. 

The goal of this work is to design the transmitter that maximizes the mutual information along N channel uses 
by jointly considering the nature of the energy harvesting process at the transmitter and arbitrary distributions of 
the input symbols, which, to the best of our knowledge, has not been yet considered in the literature. Hence, the 
main contributions of this paper are: (i.) Proving that, at the n-th channel use, the left singular vectors of the 
n-th precoder matrix are equal to the eigenvectors of the ?i-th channel Gram matrix, (ii.) Deriving an expression 
that relates the singular values of the n-th precoder matrix with the energy harvesting profile through the MMSE 
matrix, (iii.) Showing that the derivation of the optimal right singular vectors is a difficult problem and proposing 
a possible research direction towards the design of a numerical algorithm that computes the optimal right singular 
vectors. The design of this numerical algorithm is out of the scope of the current paper because our focus is to gain 
insight from the closed form power allocation that is obtained after setting the right singular vectors matrix to be 
the identity matrix and, in this scenario, the contributions are: (iv.) Deriving the optimal offline power allocation, 
named the Mercury Water-Flowing solution, and providing an intuitive graphical interpretation, which follows 
from demonstrating that the mercury level is monotonically increasing with the water level, (v.) Proposing two 
different algorithms to compute the Mercury Water-Flowing solution, proving their optimality, and carrying out an 
exhaustive study of their computational complexity, (vi.) Implementing an online algorithm, which does not require 
future knowledge of neither the channel state nor the energy arrivals, that computes a power allocation that performs 
close to the offline optimal Mercury Water-Flowing solution. 

The remainder of the paper is structured as follows. Section II presents the system model. In Section III, the 
aforementioned problem is formally formulated and solved. The graphical interpretation of the Mercury Water- 
Flowing solution is given in Section IV. The offline and online algorithms are introduced in Sections V and VI, 
respectively. In Section VII, the performance of our solution is compared with different suboptimal strategies and the 
computational complexity of the algorithms is experimentally evaluated. Finally, the paper is concluded in Section 
VIII. 
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II. System model 

We consider a point-to-point communication through a discrete-time linear vector Gaussian channel where the 
transmitter is equipped with energy harvesters. A total of TV channel uses are considered where at each channel 
use the symbol s„ G 3fJ L is transmitted. 3 

We consider that the symbols {s n }^ =1 have independent components with unit power, i.e., R s = E{s„s^} = II 
and that they are independent and identically distributed (i.i.d.) along channel uses according to Pg(s n ). As shown 
in Fig. 1, the symbol s n is linearly processed at the transmitter by the precoder matrix B„ G $t NtXL . We consider 
a slow-fading channel where the coherence time of the channel Tc is much larger than the symbol duration T s , i.e., 
T s <^Tq- Thus, a constant channel matrix H„ G $t Nr * Nt j s considered at the n-th channel use. Let K denote the 
rank of the channel matrix, i.e., K = rank(H„) = miri{TV t , TV r }, then we have that L < K. A Thus, the received 
signal at the n-th channel use is 

y n = H„B„s„ + w„, (7) 

where w„ represents the zero-mean Gaussian noise with identity co variance matrix R Wn = I^. 5 Let E„ denote 
the n-th channel use MMSE matrix, which is defined as E„ = E {(s„ — s„)(s„ — s„) T } and §„ = E {s„|y„} is 
the conditional mean estimator. 

Let us express the channel matrix as H n = Vh,^^!!^ , where A n G 5R ixi is a diagonal matrix that 
contains the L largest eigenvalues of H n and Vh„ G $l NrXL and Uh„ G $t NtXL are semi-unitary matrices that 
contain the row and column associated eigenvectors, respectively. The precoder matrix B„ can be expressed as 
B n = UB„S n V T 5n , where Ub„ G $l NtXL , £„ G 3? ixi is a diagonal matrix whose entries are given by the vector 
(T„ = [(Ti : „, . . . , <JL,n] T and Vb„ G 3? ixi is a unitary matrix. Full Channel State Information (CSI) is assumed at 
the transmitter. 

The energy harvesting process at the transmitter is characterized by a packetized model, i.e., the node is able 
to collect a packet of energy containing Ej Joules at the beginning of the ej channel access. Let J be the total 
number of packets harvested during the TV channel uses. The initial battery of the node is modeled as the first 
harvested packet E\ at e± = 1. We assume that the mean time between energy arrivals, T e , is considerably larger 
than the symbol duration time, i.e., T e 3> T s and thus we can consider that packet arrival times are aligned at the 

3 The real field has been considered for the sake of simplicity. The extension to the complex case is feasible but requires the definition of 
the complex derivative, the generalization of the chain rale, and cumbersome mathematical derivations, which is out of the scope of this work. 
Nevertheless, the extension to the complex case can be done similarly as [17] generalized the results obtained in [18]. 

4 We have considered that H n is not rank deficient, Vn, which is a realistic assumption due to random nature of the channel. 

5 Note that if the noise is colored and its covariance matrix R Wrz is known, we can consider the whitened received signal R m „ 2 y»- 
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Fig. 1. The discrete-time linear vector Gaussian channel at the n-th channel use. 
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Fig. 2. Temporal representation of energy arrivals. 



beginning of a channel use. 6 First, in Sections III-V, we consider the offline approach as it provides analytical and 
intuitive expressions. Afterwards, in Section VI, we develop an online transmission strategy where the transmitter 
only has causal knowledge of the energy harvesting process, i.e., about the past and present energy arrivals. We 
use the term pool, Tj, j = 1 . . . J, to denote the set of channel accesses between two consecutive energy arrivals. 
As in [3] and [8], we assume an infinite capacity battery since, in general, the battery size is large enough so that 
the difference between the accumulated harvested energy and the accumulated expended energy is always smaller 
than the battery capacity. A temporal representation is given in Fig. 2. 

III. Throughput maximization problem 
In this section, we study the set of linear precoding matrices {B„}" =1 that maximizes the input-output mutual 
information along N independent channel accesses, X)j=i Timer- ^( s «! Yn), where /(s„; y„) is the n-th channel use 
mutual information. The design of {B„}^ =1 is constrained to satisfy instantaneous ECCs, which impose that energy 

6 In our model, the transmitter can only change its transmission strategy in a channel access basis. Accordingly, if an energy packet arrives in 
the middle of a channel access, we can assume that the packet becomes available for the transmitter at the beginning of the following channel 
access. 
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cannot be used before it has been harvested, T s X^=i Ylner- l|BnS n || 2 < Sj=i ^j, £ = 1 . . . J. However, since in 
each pool there are several channel accesses with the same channel gains (because T c ^> T s and T e ^> T s ), instead 
of imposing the instantaneous ECCs, we can consider the mean ECCs that become T s Yfj=i 12 n er- Tr(BnBj) < 
S 3 =i Ej' ^ = 1 • • • J' which do not require prior knowledge of the transmitted symbols at each channel use as 
only the expectation of the symbols is needed. 7 

Therefore, the mutual information maximization is mathematically expressed as 

J 

max VV/(s„;y„) (8a) 

\ "/„ = i j—in^Tj 

s.t. ^ Tr(B n Bj) 1=1... J. (8b) 

j—l nerj j=l 

Before addressing the problem in (8), let us summarize the state of the art on the precoding strategy that maximizes 
the mutual information for non-harvesting nodes, which was studied in [17]-[20] and references therein. 8 In [19], 
it was shown that, in general, the mutual information, I(s n ; y n ), is not a concave function of the precoder and that 
depends on the precoder only through the matrix Z n = BjH^H„B»i- The authors of [19] also showed that the 
left singular vectors of the precoder can be chosen to be equal to the eigenvectors of the channel Gram matrix, 
i.e., Ub„ = Uh„. From this, Z n = Vr^S^A^Vg and the mutual information depends on the precoder only 
through the right eigenvectors and the associated singular values. In [18], it was shown that /(s„;y„) is a concave 
function of the squared singular values of the precoder, diag(S^), when a diagonal channel matrix is considered. 
Finally, the authors of [19] stated that the complexity in the design of the globally optimal precoder lies in the 
right singular vectors of the precoder, Vb„. Then, in [20], it was shown that I(s n ; y„) is a concave function of the 
matrix Z„ and a gradient algorithm over Z„ was derived to find a locally optimal precoder. References [18]-[20] 
considered a real channel model. The extension to the complex case was done in [17], where the authors pointed 
out that by allowing the precoder and the channel matrix to be in the complex field the mutual information can 
be further improved. Then, they proposed an iterative algorithm that determines the globally optimal precoder that 
imposes that the power constraint must be met with equality. 

When energy harvesting is considered, instead of having a single power constraint, we have a set of J ECCs 
as in (8b) and it is not straightforward to determine which of the constraints must be met with equality. This fact 

7 In general, the energy harvesting and the channel state are two independent random processes, thus, there may be situations in which only a 
few a few channel accesses separate an energy arrival from a change of the channel realization, however, note that these situations are unlikely 
since T c ^s> T s and T e 3> T s . In these improbable situations, the temporal averaging is not sufficient to ensure that the fulfillment of the mean 
ECCs implies a fulfillment of the instantaneous ECCs, however, the averaging through the different channel dimensions brings closer the mean 
and instantaneous ECCs. Thus, the mean ECCs can be used instead of the instantaneous ECCs since the cases in which they differ are indeed 
very unlikely. 

8 When there is no energy harvesting in the transmitter, the mutual information maximization problem is the one obtained after setting J = 1 
and N = 1 in (8). Thus, the mutual information is maximized for a single channel use under a power constraint. 
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implies that the algorithm introduced in [17] is no longer optimal when energy harvesting is considered. Altogether, 
(8) is not a convex optimization problem since the mutual information is not a concave function of the precoder 
and, hence, its solution is not straightforward. In the following lemma, we generalize Proposition 1 in [19] for the 
case of considering energy harvesting in the transmitter. 

Lemma 1. The left singular vectors of the n-th precoder matrix, Ub„. are equal to the eigenvectors of the channel 
Gram matrix Uh„, Vn. 

Proof: See Appendix A. ■ 
Thanks to Lemma 1, the optimal precoding matrix is B* = Uh„ S*Vg , Vn, and the dependence of /(s n ;y n ) 
on the precoder is only through E n and Vb b . In the following lines, we maximize the mutual information with 
respect to £„ for a given Vb„. By applying Lemma 1 in (7), the next equivalent signal model is obtained 

y„ = H„s„ + w„, (9) 

where H„ = Vh„ A„S„Vg and Vg is deterministic and known. To fully exploit the diversity of the channel, 
we assign the dimension of the input vector to be equal to the number of channel eigenmodes, i.e., L — K. It 
is easy to verify that the maximization of the mutual information w.r.t. S„ is not a convex optimization problem. 
However, if instead we maximize the mutual information w.r.t. the squared singular values of the precoder p„ = 
\a\ n , . . . , a\ n ] T , the obtained problem is convex, as shown in the following lines. Thus, the problem reduces to 

./ 

max ^ 5^ J(s n ; H„s„ + w n ) (10a) 

lF " ,n=1 j = l nETj 

s - t r s ^^iT Pn <^£ 3 , i = \...j. dob) 

j = l nerj j = l 

Observe that, at the n-th channel access, the input-output mutual information I(s„ ; H n s„ +w n ) is concave w.r.t. p„, 
which was proved in [18]. Therefore, the objective function is concave as the sum of concave functions is concave 
[21]. Finally, as the constraints are affine in p„, (10) is a convex optimization problem and the KKT are sufficient 
and necessary optimality conditions. In particular, the optimal solution must satisfy D Pn C = (the reader who is not 
familiar with this notation, which is presented in [1 1], is referred to [18, Appendix B] for a concise summary), where 

C is the Lagrangian that is C = liner, 7 ( s «; H«s„ + w„) - Yu=\ Pi ( T « £*=i Ener, ^Pn - Ej=i E j) . 

where are the Lagrange multipliers associated with the inequality constraints. We want to remark that in 

all the expressions derived in the remainder of the paper, n refers to some channel access contained in Tj, which 
follows from the formulation of C. In order to obtain D Pn C, we first need to determine the Jacobian matrix of the 
mutual information w.r.t. p n , which is done in the following lemma: 



Lemma 2. The Jacobian matrix of the mutual information w.r.t. p n is D Pn I(s n ; H„s n +w„) = idiag T (A^V B E„Vb„). 



Proof: See Appendix B. 
With this result, we can proceed to solve the KKT condition D Pii £ = 0: 



D p „£ = idiag T (A 2 n VlE n Vu n ) - T a £>l£ = 



=> [A* V B „E„V B „] kk = k = 1 . . . K, n e T j, (11) 

where Wj is the j-th pool water level, i.e., 

W J = ^7 ■ (!2) 

2T s 2^i =j Pi 

From (11), at each channel use, we obtain a set of K conditions that relate the power allocation in each stream 
(through the MMSE matrix) with the energy harvesting profile (through the pool's water level). Some properties 
of the water level Wj can be derived from the KKT optimality conditions: 

Pi > 0, W, (13) 

i i 

Pe( T ° - E E i) = °> V1 ( 14 > 

Plugging (13) in (12), it is straightforward to obtain the following property: 

Property 1. The water level is non-decreasing in time. 9 

From (14), we can get more insights in the solution. There are two possibilities to fulfill (14): 

• Empty Battery: This situation occurs when, at the end of the £-th pool, the node has consumed all the energy, 

i-e.. T s EU Zner, 1£p« " E?=i ^ = 0- 
> Energy Flow: This situation occurs when, at the end of the ^-th pool, the node has some remaining energy in 
the battery, which will be used in the following pools. When this happens pi = and, hence, Wi+i = We. 

Property 2. Changes on the water level are only produced when at the end of the previous pool the node has 
consumed all the available energy. 

Note that the ECCs take into account the energy spent by the node over all the dimensions. Thus, these two 
properties also hold in a scalar channel model as proved in [6, Theorem 3]. 

Since the problem in (10) is convex, by using (11) and Properties 1 and 2, we can construct efficient numerical 
algorithms to compute the optimal power allocation, {p*}^ =1 , for a given Vb„- The maximization of the mutual 

'This property is only valid under an infinity battery capacity assumption. When a finite battery is considered the water level may increase 
or decrease [7]. 
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information w.r.t. Vb„ is indeed much more complicated as pointed out in [19] for the non-harvesting scenario. 
In this context, in this work, we focus on the particular case in which Vb„ = Ik because, in spite of not being 
necessarily the globally optimal precoder, it leads to an analytical closed form power allocation that allows an 
intuitive graphical representation of the solution, as it is explained in the next section. Observe that for any other 
choice Vb„ ^ lit, we must resort to numerical methods to compute the optimal power allocation. 

The design and development of a numerical algorithm that computed the globally optimal precoder at each 
channel access would be an interesting research problem in its own and is left for future research. We believe that a 
possible starting point would be to analyze how to expand the algorithms presented in [20] and [17], which exploit 
the concavity of the mutual information w.r.t. the matrix Z„ and the fact that the power constraint must be met with 
equality, to the energy harvesting scenario. Note that if we knew the optimal total power allocation in each channel 
access, we could run N times the algorithm proposed in [20] to obtain the globally optimal precoder in each channel 
access, however, this approach has two major drawbacks. First, the optimal total power allocation in each channel 
access is not known a priori and its computation is not straightforward since the total power consumptions of the 
different channel accesses belonging to the same pool must simultaneously satisfy (11). The second drawback is 
the required computational burden since any iterative approach requires a new estimation of the MMSE matrix, 
E„, at every iteration since it depends on Vb„ and £„. These two reasons makes challenging the applicability 
of the proposed approach and, hence, different alternatives to find the globally optimal precoder may be required. 
Altogether, we believe that the development of a numerical algorithm that computes the globally optimal precoder 
for a WEHN is the object of a new paper in its own and is left for future research. 

IV. The Mercury Water-Flowing solution 

In the remainder of the paper, we consider a communication system in which the precoder is constrained to satisfy 
Vb„ = I if or, equivalently, a communication system such that both the precoder and channel matrices are diagonal. 
In spite of the fact that total achievable mutual information is reduced by forcing Vb„ = Ik, we consider that it 
is interesting to study this scenario for the following three reasons: (i.) The system y' n = A n S„s„ + w^, with 
being the observed noise at the receiver, is commonly encountered in practical systems where, for simplicity at the 
decoder, independent symbols are transmitted in each dimension (e.g., in multi-tone transmissions like Orthogonal 
Frequency Division Multiplexing (OFDM)), and it has been broadly considered in the literature, indeed, the "H g WF 
solution was derived for such an input-output system model in [14]. (ii.) The optimal power allocation, which is 
named Mercury Water-Flowing, accepts a closed form expression and an intuitive graphical representation, (iii.) 
We believe that the intuition gained thanks to the Mercury Water-Flowing graphical interpretation may help for the 
design of the algorithm that computes the globally optimal precoder of the problem in (8). 
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In this context, the input-output model y' n = A„S„s n + w' t can be obtained from the general model in (9) by 
setting Vg = Ijj and y' n = Vjj y„. From this, we obtain that the equivalent noise is w' n = Vjj w„. Thus, 
a set of K independent parallel streams are observed at each channel use. The received signal in the fc-th stream 



i s v'k n = ^/^fe,n cr fc n s k,n + w' k n , where the transmitted symbol is the k-th component of s„, i.e., Sk. n = [s n ]k, 
w 'k n = [ w 'iJfc i s tne observed noise, Xk.n = [^n]kk is the channel gain, and cri n = [E^]kfc is the transmission 
radiated power. Therefore, in this section we solve the following optimization problem: 

J 

max y y I(s n ;y' n ), subject to (10b). (15) 
at. „ . ~r 

Note that /(s n ;y^) = /(s„;y„) since a linear unitary rotation in the received signal does not affect the input- 
output mutual information [17]. Thus, the power allocation that maximizes (15) is equal to the one that maximizes 
(10) and it can be obtained by particularizing (11) with Vg n = Ik, i.e., [E„] fcfc = Afc 1 ty . . From where, it follows 
that 



1 



k,n A 



-mrase. 



k. 



{ min { h w^})' Vfc ' V ^' V " G ^ (16) 



where mmse^ (■) is the inverse MMSE function, defined as in [14], that returns the SNR of the fc-th stream for 
a given MMSE, which depends on the probability distribution of Sk, n - 

To present the graphical interpretation of the solution, we need to reformulate (16) as 

1 „ ( 1 



where Gk(ip) is defined in (6) [14], depends on the modulation used, and satisfies the next lemma: 
Lemma 3. The function Gfc(?/>) is monotonically decreasing in tp. 

Proof: See Appendix C. ■ 

Remark 1. To demonstrate the validity of the graphical representation presented in this section, we need to 
analytically demonstrate that Gfc('0) is monotonically decreasing in ip. In [14], it was already stated that Gk(-) is 
decreasing, however, the authors did not provide an analytical proof for their statement. Therefore, we consider that 
Lemma 3 and its explicit proof are crucial to validate the graphical representation introduced in this section. 

Observe the similarity of the power allocation found in (17) with the H g WF in (5) [14]. The main difference 
of our solution is that, due to the nature of the energy harvesting process, the water level depends on the channel 
access. Indeed, from Properties 1 and 2 we have seen that the node is able to increase the water level as energy is 
being harvested. 

Moreover, observe that if we particularize (17) for Gaussian distributed inputs, which have Gk(ip) = 1, V-0, (see 
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[14]), the DWF solution in [7] is recovered. Therefore, the mercury factor gives a measure of how power allocation 
is modified when using non-Gaussian input distributions. 

Let - Hg fe,n ^(Wj) be the mercury level of the k-th stream at the n-th channel use, i.e., 



which depends on the gain and water level of the channel use. Then, the power allocated in a certain stream is the 
difference between the water and mercury levels, i.e., a\ n = (]¥j — 'Hg k ' n ^ (Wj)\ . The solution interpretation 
presented in this section is based on the fact that the mercury level is monotonically increasing in Wj, which 
follows directly from Lemma 3, and generalizes both the "H g WF and the DWF solutions derived in [14] and [7], 
respectively. The Mercury Water-Flowing interpretation, depicted in Fig. 3, is the following: 

1) Each parallel channel is represented with a unit-base water-porous mercury-nonporous vessel 10 . 

2) Then, each vessel is filled with a solid substance up to a height equal to A^. 

3) A water right-permeable material is used to separate the different pools. 

4) Each vessel has a faucet that controls the rhythm at which mercury is poured. The faucet modifies the mercury 
flow so that the relation between mercury and water levels in (18) is always satisfied. 

5) Simultaneously, 

• The water level is progressively increased to all pools at the same time, adding the necessary amount of water 
to each pool. The maximum amount of water that can be externally added at some pool is given by the pool's 
harvested power (Ej/T s ). Let the water freely flow right through the different pools. 

• Mercury is added to each of the vessels at a different rhythm which is controlled by the vessel's faucet. 

6) The optimal power allocation in each parallel channel is found when all the pools have used all the harvested 
energy and is obtained as the difference between the water and mercury levels. 



We have designed two different algorithms to compute the optimal Mercury Water-Flowing solution, namely, the 
Non Decreasing water level Algorithm (NDA) and the Forward Search Algorithm (FSA), which are presented in 
Sections V-A and V-B, respectively. Afterwards, in Section V-C, we prove the algorithms' optimality and analyze 
their computational complexity. 

As shown by the KKT optimality conditions, the water in a certain pool may flow to the pools at its right (i.e., 
from prior to later time instants). This way, the water level over a consecutive set of pools may be equalized. This 
set of constant water level pools is referred to as an epoch, £ m ,m = 1 . . .M, where M is the total number of 

'"The vessel boundaries are not depicted in Fig. 3 for the sake of simplicity. 




(18) 



V. Mercury Water-Flowing offline algorithms 
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epochs and it is unknown a priori. Note that, since the pools are a partition of the epochs, a certain pool Tj is only 
contained in one epoch. However, an epoch may contain several pools, therefore, M < J. 

To compute the power allocation in (17), we just need to determine which pools are contained in each epoch 
£ m as, once the epochs are known, the optimal power allocation of the epoch can be found by performing the 
Mercury AVaterfilling Algorithm ("HgWFA) introduced in [14], where the m-th epoch water level, W m , is found by 
forcing that the energy expended in the epoch has to be equal to the energy harvested, which follows from Property 
2. 

The following two algorithms use a different approach two determine the epochs: 
A. NDA 

The NDA uses the fact that a water level decrease is suboptimal, which follows from the KKT conditions (see 
Property 1), to compute the optimal power allocation as follows: 

1) Initially, set M := J, i.e., every epoch contains one pool £ m := {r m }, m = 1 . . . M. 

2) Perform the HgWFA in [14] to every epoch to obtain the water level, W m> in each epoch. 

3) Look for some epoch, m', at which the water level decreases, i.e., W m > > W m >+i'. 
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• If some epoch is found, merge this epoch with the following epoch, i.e., £ m > := £ m > U £ m '+i- The harvested 
energy of the resulting epoch is the sum of the two original epochs. Then, the total number of epochs has been 
reduced by one, i.e., M := M — 1. Perform the H g WFA to obtain the new water level of the m'-th epoch, i.e., 
W m i, and go back to 3. 

• If no epoch is found, the optimal M has been found along with the optimal power allocation. 

B. FSA 

The FSA determines the different epochs by finding the optimal transition pools, {7^}„ =1 , that are defined as 
the first pool of each epoch. As stated before, once the epochs are known the optimal power allocation is determined 
by applying the 'H g WFA to each epoch. 11 To determine {7^}^ =1 , we have designed a forward-search algorithm 
that extends the algorithm introduced in [6] to take into account arbitrary input distributions. We explain how to 
obtain 7^* and the others are found in the same manner: 

1) Assume that the first epoch contains all the pools, S\ := {ri,T2 ■ • • ,tj}. 

2) Perform the "HgWFA in [14] to the epoch. 

3) Check whether all the ECCs within the epoch are fulfilled: 

• If they are not fulfilled, remove the last pool from E\ and go back to step 2. 

• If they are fulfilled, the optimal transition pool, 7^*, is the first pool not included in the epoch. 

The same procedure is repeated to determine the following transition pools until the A-th pool is included in 
some epoch. When this happens, the optimal power allocation has been found for all the channel accesses and 
streams. 

C. Optimality and performance characterization of the offline algorithms 

In this section, first, we demonstrate the optimality of the NDA and the FSA, which is presented in Theorem 1 
and, afterwards, we characterize their associated computational complexity. 

Theorem 1. Both the NDA and the FSA compute the optimal power allocation given in (17). 

Proof: See Appendix D. ■ 
With the previous theorem, we have demonstrated that both algorithms compute the optimal power allocation, 
however, the computational cost of such a computation may be very different. To evaluate this, in Appendix E, we 
have conducted an exhaustive study on the computational complexity of each of these two algorithms. 

"Observe that, by definition, 7j* is the first channel use. 
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Table I 

Computational complexity of the NDA and the FSA in the best, worst and average case scenarios. 



Our performance analysis is three-fold, namely, the best, worst, and average computational complexities are 
computed. Note that both algorithms internally call the H g WFA a certain number of times to find the optimal 
solution. Let CV^wfa denote the number of calls to the W S WFA required to compute the Mercury Water-Flowing 
solution, which depends on the algorithm itself and on the dynamics of the energy harvesting process. In this context, 
the best or worst computational complexity is the performance when the minimum or maximum number of calls to 
H g WFA are required, respectively. The average computational complexity uses a probabilistic model to compute 
the average number of calls to 'HgWFA. Basically, for the NDA we assume that there is a fixed probability q that 
the water level decreases from epoch to epoch, whereas, for the FSA we assume that there is a fixed probability p 
that a certain ECC is not satisfied. Both p and q can be experimentally adjusted depending on the energy harvesting 
profile. The computational complexity in terms of operations (Op.), as well as, in terms of C-h s wfa is summarized 
in Table I, where a is a constant parameter that depends, among others, on the size of the MMSE table required 
to compute the inverse MMSE function and on the tolerance used in the stopping criteria of the % g WFA. The 
details of the derivations of the different computational complexities can be found in Appendix E. In Section VII-B, 
the theoretical results on the algorithms' computational complexities are compared with the ones obtained through 
simulation. 

VI. Online algorithm 

Up to now, we have assumed that the transmitter has non-causal knowledge of both the CSI and the energy 
harvesting process, which is not a realistic assumption in practice. Therefore, the Mercury Water-Flowing solution 
provides an upper bound on the achievable mutual information of practical schemes in which Vb„ = Ik- In 
this section, we develop an online algorithm, which is strongly based on the optimal offline solution, the Mercury 
Water-Flowing power allocation, but that does not require future knowledge of neither the energy arrivals nor the 
channel state, that computes a suboptimal power allocation of the problem in (15). 

Let F w be the flowing window that is an input parameter of the online algorithm that refers to the number of 
channel accesses in which the water is allowed to flow, which can be obtained by a previous training under the 
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considered energy harvesting profile, and let an event denote a channel access in which a change in the channel 
state is produced or an energy packet is harvested, i.e., s t = {n|A n _i ^ A n } U {n\n = e,}, t = 1 . . .T, where 
T € [J, N]. In this context, the proposed online algorithm proceeds as follows: (1.) The initial energy in the battery, 
Ei, is allocated to the different streams of the first F w channel accesses according to the H g WFA where the 
channel is expected to be static and equal to the gain of the first channel use A n = Ai,Vn € [1,^]. (2.) When 
the transmitter detects an event, it updates the allocated power of the channel accesses n e [,s t , min{ s t + F w — 1 , N}] 
by using the H g WFA with the remaining energy in the battery and with the energy of the harvested packet (if 
the event is an energy arrival), i.e., Y^j\ St <e - ~ T s J2n=i Sfe a \ n > an< ^ by assuming that the channel remains 
constant during the flowing window, i.e., A„ = A St , Vn £ [st,vnm{st + F w — 1,N}]. Note that the transmitter 
may stay silent in some channel accesses if the difference between two consecutive incoming energy packets is 
greater than the flowing window, ej — ej-i > F w . 12 Step (2.) is repeated until the iV-th channel access is reached. 
The proposed online algorithm satisfies ECCs and, as pointed out, does not require future information of neither 
the channel state nor the energy arrivals. 

The performance in terms of achieved mutual information depends on the correctness of the estimation of the 
flowing window, F w , as discussed with the numerical analysis in Section VII. In summary, this online algorithm 
provides us a lower bound on the mutual information that can be achieved with sophisticated online algorithms that 
make use of precise statistical models of the energy harvesting process and channel state. 13 

VII. Results 

This section first evaluates the gain of the proposed Mercury Water-Flowing solution with respect to other 
suboptimal solutions and, secondly, it presents an analysis through simulation of the computational complexity of 
the NDA and FSA. 

A. Results on Mercury Water-Flowing solution 

In this section, we evaluate the mutual informations obtained with the optimal offline solution, the Mercury 
Water-Flowing CH g -WFlow), and with the online policy presented in Section VI. 

To the best of our knowledge, there are no offline algorithms in the literature that maximize the mutual information 
by jointly considering energy harvesting at the transmitter and arbitrary distributed input symbols. In this context, 
we use the following three algorithms, which are optimal in different setups and have been adapted to the energy 

12 This situation rarely takes place in practice since, in most common situations, F w is several times the mean number of channel accesses per 
pool. For example, in the simulated framework presented in Section VII, we have obtained that F w is 4.4 times the mean number of channel 
accesses per pool. 

13 A myriad of works have dealt with channel modeling, however, having a precise statistical model of the energy harvesting process is indeed 
not trivial as it depends on many factors such as the harvester used by the node (e.g., a solar panel, piezoelectric generator, etc.), the node's 
placement, mobility, etc. 
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harvesting scenario, as a reference to evaluate the mutual information achieved by the proposed offline and online 
solutions: (i.) The DWF solution in (1) that is the optimal offline power allocation for a WEHN when the distribution 
of the input symbols is Gaussian, (ii.) Pool-by-Pool Waterfilling (PbP-WF) that uses the WF power allocation in 
(2) by forcing that the harvested energy in a certain pool is expended in the channel accesses of that same pool, 
(iii.) Pool-by-Pool Mercury /Waterfilling (PbP-H g WF) where the power allocation is obtained by using the H g WF 
solution in (5) and forcing that the harvested energy in a certain pool is expended in the channel accesses of that 
same pool. 

We have considered a channel matrix of rank K = 4, where the channel gains are generated randomly. The 
modulations used in each stream are BPSK, 4-PAM, 16-PAM, and 32-PAM, respectively. The symbol duration is 
T s = 10 ms and N = 100 channel accesses have been considered during which a total of J = 40 energy packets 
are harvested. Energy arrivals are uniformly distributed along the channel accesses and with random amounts of 
energy, which are normalized according to the total harvested energy that varies along the a>axis of Fig. 4. The 
y-axis shows the mutual information obtained with the different strategies. After some training in this scenario, we 
have obtained that the optimal flowing window is F w = 11 channel accesses. As shown in Fig. 4, our proposed 
solution, the 'Hg-WFlow, outperforms all the suboptimal strategies. The improvement of the Hg-WFlow w.r.t. the 
PbP-HgWF comes from letting the water to flow across pools and, hence, it directly depends on the parameter 
J since the higher is the number of pools, the higher is the mutual information gain that can be achieved by 
letting the water flow. 14 The same happens with the improvement of the DWF w.r.t. PbP-WF. On the other hand, 
the mutual information gain of the PbP-HgWF and Hg-WFlow w.r.t. their respective WF strategies, PbP-WF 
and DWF, comes from the use of mercury in the resource allocation. Thus, when the energy availability is low, 
both perform similarly because the node is working in the low SNR regime in which the mutual information of 
finite alphabets is well approximated by the mutual information of the Gaussian distribution [22]. However, when 
the energy availability is high, the PbP-HgWF and "H g -WFlow achieve a higher mutual information than their 
respective WF strategies since the mutual information of finite constellations asymptotically saturates (not more than 
log 2 Q bits of information can be sent per channel use). Finally, note that, in spite of not having knowledge of the 
energy arrivals nor channel state, the online power allocation performs close to the the offline optimal "H g -WFlow 
in the low SNR regime. When the available energy increases, the gap between the Mercury Water-Flowing and the 
proposed online algorithm also increases, nevertheless the online algorithm still presents a reasonably good mutual 
information outperforming any Pool-by-Pool strategy. 

The study of the performance in the static scenario is of special interest because the assumption of having future 
knowledge of the channel state, which has been used for the design of the optimal offline solution, becomes realistic 

14 When J = 1, the solid and dashed curves overlap since there is only one pool. 
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Fig. 4. Mutual information for the different transmission strategies versus total harvested energy. 



when the channel is static. We have evaluated the achieved mutual information in the above setup for the static 
channel case and we have obtained similar results than the ones in Fig. 4, where the only difference is that the 
achieved mutual information of the different algorithms in the static case is slightly lower since there is less channel 
gains diversity to assign the available energy. 15 

In Fig. 5, the power allocation obtained by the Mercury Water-Flowing solution in a single simulation is shown 
for N = 20 and K = 4, where the modulations used in the streams 1-4 are BPSK, 4-PAM, 16-PAM and 32-PAM, 
respectively. Six energy arrivals are produced at the beginning of the channel accesses marked with a triangle. The 
gains have been generated randomly along channel uses, but fixed constant along streams to ease the observation of 
the mercury level obtained for the different modulations. As expected from Property 1, the obtained water level is 
an increasing step-wise function. Observe that the solution contains three epochs, i.e., three different water levels, 
where the pools contained in each epoch are E\ = {ti}, £2 = {T2, T3, T4}, and £3 = {T5, tq}. Moreover, observe that 
under the same channel gain and water level, the mercury level decreases as the modulation dimension increases. 

15 The figure of the static scenario has been omitted for the sake of brevity. 
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Fig. 5. Graphical representation of the Mercury Water-Flowing solution. The red, gray, and blue solid bars represent the inverse of the channel 
gain, the mercury and the water levels, respectively. The allocated power is obtained as the difference between the water level and the mercury 
level. 



B. Results on the algorithms ' performance 

In Section V-C, we have given a summary of the computational complexity of the NDA and FSA (Table I 
summarizes the obtained results). In this section, we compare the theoretical and experimental performance of both 
algorithms. 

From the simulations, we confirm that, in the best and worst case scenarios, the experimental computational 
complexity shown in Fig. 6 fits the theoretical results presented in Table I. Regarding the average case scenario, 
the mean number of calls to H g WFA of the NDA fits the analytical expression E {C'^ i ^ 1 FA } = J(q + 1) — q 
for a value of q = 0.98. Regarding the FSA, the mean obtained through simulation and the analytically computed 
expression E{C^ FA } = (£ + f - l) p+1 differ from one another. Observe that the quadratic and linear terms 
of J have the same weight independently of the value of p. However, it is easy to observe in Fig. 6 that the linear 
component dominates over the quadratic. Therefore, there is a mismatch between the analytical and experimentally 
obtained expressions. We believe that this mismatch is due to the fact that in order to obtain some tractable model 
(see Appendix E), we have assumed that all the ECCs have the same probability p of not being satisfied, however, 
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in reality this probability is not necessarily equal but depends on the dynamics of the energy harvesting process. 

Regardless of the aforementioned mismatch, we observe that, in our simulated energy harvesting set up (the 
amount of energy in the packets is uniformly distributed), both algorithms have a similar performance in the 
average case scenario. Note that the difference between the best and worst case scenario is much smaller for 
the NDA than for the FSA. This comes from the fact that, in the worst case scenario, the FSA has a quadratic 
dependence in J, whereas, for the NDA the dependence is linear. This makes the NDA more robust in front of 
changes in the energy harvesting profile. In other words, if the energy harvesting profile changes, the FSA has more 
margin to either improve or degrade its performance. For instance, if the node initial battery is very high and the 
node is operating in the sunset (the amount of harvested energy at the beginning of the transmission duration is 
higher than the amount harvested at the end) it is likely that the performance of the FSA is close to the best case 
scenario, i.e., a single call to the "H g WFA. On the other hand, if the battery is almost empty at the beginning and 
the node operates in the sunrise the performance of the FSA will be very poor. 

To conclude the discussion between the NDA and the FSA, we want to highlight again that the NDA is more 
robust to changes in the the energy profile characteristics. However, the FSA may be preferable in certain energy 
harvesting profiles as in its best case performace just requires a call to the H S WFA. Therefore, we believe that 
the algorithm selection must be done by taking into account the energy harvesting profile and the environmental 
conditions in which the node is operating. 

VIII. Conclusions 

In this paper, we have considered a WEHN transmitting arbitrarily distributed symbols through a discrete-time 
linear vector Gaussian channel. We have studied the precoding strategy that maximizes the mutual information 
by taking into account causality constraints on the use of energy. We have proved that the optimal left singular 
vectors of the precoder matrix diagonalize the channel, similarly as in the optimal precoder for the case of non- 
harvesting nodes. We have derived the expression [A^V^EnVsJ^ = ^ that relates the singular values of 
the precoder (through the MMSE matrix) with the energy harvesting profile (through the different water levels). 
The derivation of the optimal right singular vectors, Vg , is left as an open problem. Then, we have derived 
the Mercury Water-Flowing solution, the optimal power allocation when Vb„ = Ik, which can be expressed 
in closed form and accepts an intuitive graphical interpretation based on the fact that the power allocation in a 
certain stream is the difference between the water level and the mercury level, which, as shown in this paper, is 
a monotonically increasing function of the water level. Additionally, we have developed two different algorithms 
that compute the Mercury Water-Flowing solution and we have analytically and experimentally evaluated their 
computational complexity. We have also proposed an online algorithm that only requires causal knowledge of the 
energy harvesting process and channel state. Finally, through numerical simulations, we have shown a substantial 



21 



Performance of the algorithms 




J (N = 300) 
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increase in the mutual information w.r.t. other suboptimal offline strategies, which do not account for the shape, size 
and distribution of the input symbol or do not exploit the water level equalization across pools, and we have seen 
that the mutual information achieved with the online algorithm is close to the one of the Mercury Water-Flowing 
solution. 

Appendix 

A. Proof of Lemma 1 

Let us assume that the optimal precoding matrices of the channel accesses n = 2 . . . N are known, i.e., {B* }^_ 2 . 
Then, we focus on finding the optimal precoding matrix of the first channel use B^. The problem in (8) is equivalent 
to 

max J(si;yi) + a (19a) 

Bi 

I 

s.t. T s Tr(BiB[) + b + c(£) < ^Ej, £=1...J, 

3 = 1 
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£ p „H„ = -(V B „ ® V^A^SKDiagK; 1 ) = -(V B „ ® V H „ A„)S A -Diag(diag(i:- 1 )) 

= i(V B „®VH„A„)SKSl.(lK®S- 1 )S A = i(V B „®VH„A„)(Iif®S- 1 )S /f , (21) 

D p „/(s„;H„s„+w„) = ivec T (H„E n )(V B „«)VH„A n )(l fc ®S; 1 )SK (22) 

= ivec T (H„E„)(V B „®V H „A n S- 1 )S K (23) 

= ivec T ((V H „A n S j ; 1 ) T H Il E„V B „)s K (24) 

= ivec T (A2V B „EnV B „)S A . = idiag T (A2V Bi E n V B J, (25) 



where a, b and c(£) do not depend on Bi. By only keeping the most restrictive constraint, which is denoted by 
Pc, the previous optimization problem reduces to 

max 7(si;yi) (20a) 

Bi 

s.t. Tr(BiBf) < P c - (20b) 

Finally, once the problem is expressed as (20), it is known from [19, Pip. 1] that the left singular vectors of BJ 
can be chosen to coincide with the eigenvectors of Rf^, i.e., U Bl = Uhj- A similar approach can be applied to 
show that {B*}^ =2 diagonalize their respective channels. □ 

B. Proof of Lemma 2 

By applying the chain rule, we have that D Pn /(s„; H„s n + w„) = Dpj I(s n ; H„s„ + w„) D Pn H„. The first 
term in the previous equation can be easily derived from (4) as Dg I(s n ; H„s„ +w„) = vec T (H n E n ). The second 
term, D Pii H„, is given in (21), where the first equality can be proved in a similar manner than D^P in [18, Proof 
of Theorem 5]. Sk £ ^ K xK is the reduction matrix introduced in [18] (See Appendix F for a concise summary 
on the properties of Sk)- In the third and fourth equalities, we have applied Properties 6 and 8 in Appendix F, 
respectively. Therefore, D Pri I(s n : H„s„ + w„) is derived in (22)-(25), where, in (23) and (24), we have used that 
(A ® B)(C g> D) = AC <g> BD and vec(ABC) = (C T ® A)vecB for any matrices A, B, C, and D such that 
the matrix products AC, BD, and ABC are well defined [11]. Finally, (25) follows from the definition of the 
reduction matrix (See Appendix F). This concludes the proof. □ 
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C. Proof of Lemma 3 
Let ip be some fixed MMSE that can be obtained as 

ip = mmsecisnro) — mmseA(snrA), (26) 

where mmseoisnrQ) and mmseA(snrA) give the MMSE as a function of the SNR for a Gaussian and for an 
arbitrary input distribution, respectively. Thus, snrc and snrA are the associated required SNRs to achieve the 
error ip for these distributions. 

Similarly, the required SNR to obtain a certain error can be computed by the inverse MMSE function as sure = 
mmse^}{ip) and snrA = mmse^ 1 (ip) . 

For the Gaussian case, it is broadly known that ip = mmseMsnra) = tt-^ — with derivative imm ^ e o( snr G) _ 
(l+snrc)" - Similarl y' snr G = mmse^ 1 ^) = i - 1 and dmm ^ W = 

Note that for any generic function f(x) it is verified that d ^ j/^^ = d ^ g/Tx) ^dx^ = L By applying the 
previous property, the following relation is obtained: 

dmmse^ 1 (ip) dmmsea(snrc) dmmse^ 1 (ip) dmmseA(snrA) 
dip dsnrc dip dsnrA 

Recall that G(ip) = ^ — mmse^ 1 (ip) as ip € [0, 1]. Then, its derivative is 

AG(ip) —1 dmmse^ 1 ^) dmrnse^ 1 ^) dmmse~^ 1 (ip) 

dip ip 2 dip dip dip 

di/> / dmmse a (snrA) dmmseo(snrQ) 



Ammse A (snr A ) ^ dsnrA dsnr G 

In [23], it was recently shown that mmseA(snrA) = E{M2} and dmm ( ^ r ^" r ' A - ) = —ElM^}, where M2 = 
var{x\y/ snr ax + n}. Therefore, the first term of the previous equation is always positive since both the MMSE 
and the inverse MMSE functions are decreasing for any distribution. In (27)-(30), we show that the second term 

is positive, where in (29), we have used that snrn = —, r — 1, which follows from (26). In (30), we 

have used the recently found expressions of the MMSE and its derivative [23]. Finally, as the variance is positive, 



24 



iG ffl < and G(ip) is a monotonically decreasing function. □ 
D. Proof of Theorem 1 

The optimality of the algorithms is proved by demonstrating that the power allocation obtained by means of each 
of the algorithms satisfies the KKT sufficient optimality conditions: 
(1.) ^ = 0,Vfc,n. 

(2.) t s y,U E„ erj £f=i <„ < E U Ej, t = i...J. 

(3.) p t >Q,l=l...J. 

(4.) Pi(T s EU ^ner 3 Etl <n ~ E "=1 ^) = CU = 1 . . . J. 

Moreover, we know that by the end of the transmission the battery must be empty since, otherwise, the remaining 
energy in the battery can be used to increase the total mutual information. Thus, (2.) must be met with equality for 
£ = J. Note that both algorithms compute a power allocation strategy that satisfies ECCs and that by the end of 
the last channel access all the energy has been used. Therefore, (2.) is satisfied W and it is satisfied with equality 
for 1 = J. From Property 1, if the water level is non-decreasing in time then (3.) can be verified. In the NDA, the 
water level is clearly non-decreasing in time. Regarding the FSA, if some ECC is not satisfied, it is because the 
water level must be reduced before the point where the ECC is not satisfied and increased afterwards. Indeed, this is 
what the algorithm does in the procedure of finding the optimal epochs. Therefore, (3.) is also satisfied in the FSA. 
Finally, since both algorithms compute the optimal power allocation within an epoch by using the "H g WFA, where 
the water level is found by forcing that all the available energy must be used by the end of the epoch, conditions 
(1.) and (4.) are satisfied. With this, we have demonstrated that the power allocation computed by the NDA and 
the FSA is the optimal power allocation. □ 



E. Computational complexity of the algorithms 

In this appendix, we study the performance of the two algorithms that compute the Mercury Water-Flowing 
solution, the NDA and the FSA. 

We have carried out a three-fold analysis, namely, the best, worst and average computational complexity. As 
mentioned before, both algorithms internally call the "H g WFA a certain number of times to find the optimal 
solution. The performance is evaluated in terms of operations and number of calls to the H g WFA required to 
compute the Mercury Water-Flowing solution, C-^wfa- 

Before getting into the complexity of each of the aforementioned scenarios, let us first compute the complexity of 
the 'HgWFA when the algorithm computes the power allocation of NK parallel channels, where N and K denote 
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the number of channel accesses and streams, respectively, i.e., 

CCn sWF (N,K) = aNK, (31) 

where a is a constant parameter that depends, among others, on the size of the MMSE table required to compute 
the inverse mmse function mmse'^ 1 { ) and on the tolerance used in the stopping criteria of the HgWF. Now, let 
us proceed to compute the computational complexity of the NDA and FSA. 

1 ) Computational complexity in the best case scenario: 

NDA: The best case scenario for the NDA occurs when the resulting water-levels of applying the HgWFA at 
each pool are non-decreasing throughout all the transmission. Thus, the best case computational complexity for the 
NDA is 

J J 
CC§ DA (N, K, J)=J2 CC Hs wf{L 3 ,K) = aL K = aNK, (32) 
i=i i=i 

where Lj is the number of channel accesses contained in Tj and, accordingly, ^2j=i Lj = N. Note that the number 
of calls to the T-L g WF is C-h s wfa = J- 

FSA: Regarding the FSA the best performance is obtained when the algorithm can stop at the first iteration, i.e., 
after applying the % g WFA to the N channel accesses it is observed that the resulting power allocation satisfies all 
energy causality constraints, i.e., 

CC$ SA (N, K, J) = CC Hs wf{N, K) = aNK. (33) 

Note that the number of calls to the H e WF for the FSA in the best case scenario is Cu s wfa = 1. 

Observe that, even though C-H g wFA differs from one algorithm to another one, they achieve the same computa- 
tional complexity in terms of operations in the best case scenario. However, note that the best case scenario for the 
FSA occurs when the water level of the optimal power allocation remains constant throughout all the transmission 
time, in other words, there is a single epoch. However, the best case scenario for the NDA is completely the 
opposite, the water level is different at every pool and, thus, the total number of epochs is J. 

2) Computational complexity in the worst case scenario: 

NDA: The worst case computational complexity for the NDA is produced when at every iteration of the algorithm 
it is observed that the water level is decreasing in some pool transition. Fig. 7 shows an example of how the algorithm 
proceeds for J = 4. In the first iteration a total of J calls to H g WF are required. Then, in the second iteration, an 
additional call is performed to merge the first two epochs where it is observed that the water level is decreasing. As 
we are considering the worst case scenario, the resulting water-levels will be decreasing at some epoch transition 
and an additional call is required until all pools have been merged in a single epoch, therefore, the worst case 
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Fig. 7. Representation of the NDA algorithm. 
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computational complexity for the NDA is 



CC% DA (N,K,J) 



CC u%WF {L h K) + CC Hg w F (jL 3 ,K) 

3 = 1 3=2 
J 

aKN + Y aKjN/J 

3=2 



(34) 

(35) 
(36) 



where the first summation comes from the first iteration of the algorithm and the second one comes from merging 
the pools with decreasing water level, i.e., iterations from 2 to J. In (35), we have made the simplification of having 
equal length pools, i.e., Lj = N/J, Vj. The number of calls to H g WF is Ch^wfa = 2 J — 1. 

FSA: The FSA starts by assuming that the first epoch contains all the pools, then, it performs H g WF and 
checks whether the energy causality constraints are satisfied, which are not as we are considering the worst case 
scenario. Then, it removes the last pool from £\ and tries again and so forth until E\ just contains one pool and 
then the constraints must be satisfied. Therefore, a total of J iterations are required to determine £*. Similarly, 
J — 1 iterations are required to determine The computational complexity at each iteration is summarized in 
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Iteration 


Epoch 


1 r\tnnlAvit\' 
V^UlllUlCAlLV 


1 


Ex = {n, ...,tj} 




2 


£l = {Tl,...,Tj_i} 




J 


Si = {n} 






Total £\ = a,K\JL\ + (J - l)L 2 + • • • + Lj] 



J+l 
2J-1 


£2 = {t 2 ,...,tj} 
£ 2 = {r 2 } 


aKL 2 




Total = aK[(J - 1)L 2 + ( J - 2)L 3 + ■■■ + Lj] 



J(J+l) 

2 


£j = {tj} 


aKLj 




Total £*j = aKLj 



Table II 

Computational complexity of the FS A in the worst case scenario. 



Table II from where we can conclude that the worst case computational complexity of the FSA is 

J 

CCf SA {N,K,J) = J2 aKL 3( J -i + 1 ) ( 37 > 

3=1 

= aKjY.j J -i 2+ j = °^ KNj2 ^ ( 38 > 

3=1 

where in (38) we have made the simplification of having equal length pools, i.e., Lj — N/J, Vj. As every iteration 
performs a call to H g WF, the total number of calls is Cu s wfa = ^ ■ 

3) Computational complexity in the average case scenario: For the average case scenario, due to the inherent 
difficulty of determining the computational complexity measured in operations, we have just derived the complexity 
in terms of calls to the H g WFA, i.e., C7^ g wFA- By doing this, we can see how the computational complexity is 
affected by the number of energy arrivals J. 

NDA: We start by analyzing the average performance of the NDA. Let qj, j = 1 ... J — 1, be the probability 
that the water-level decreases at some pool transition. Let us assume equal probability at all pool transition qj = q, 
Vj. Let C|)^w FA be a random variable that, for a certain call to the NDA algorithm, denotes the number of calls 
to the 'HgWFA. Note that the minimum number of calls to the 'HgWFA is J and, from here, an additional call 
is produced every time that a water level decrease is produced. Observe that this additional number of calls is a 
binomial distribution of parameters J — 1 and q, i.e., B (J — 1, q). Therefore, C^H^ FA = J + B (J — 1, q) and the 
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mean and variance are 

M^wfaI = J + E{E(J-l,q)} = J+{J-l)q=J(q+l)-q, (39) 
var{C#3& A } = var{B (J - l,q)} = (J — 1)<?(1 — ?). (40) 

FSA: Similarly for the FSA, let Pj, j = 1 ... J — 1, denote the probability that the j-th energy causality constraint 
of the FSA is not satisfied. We assume that this probability is equal for all the constraints Pj = p, Vj. Let ^ FA 
be a random variable that, for a certain call to the FSA algorithm, denotes the number of calls to the "H g WFA. To 
determine E {C^wfa} f° r a general J, let us first obtain ^fa f° r some specific values of J as a function of 
the broken constraints. Note that up to J — 1 constraints can be broken. In Tables III, IV, V, t/ and X denote that 
a certain constraint is satisfied or broken, respectively. For example, Table III shows C^ S ^ FA when J = 3 and the 
energy constraints that can be broken are in the transitions of t\ — > t-x, which is depicted in the first column, and 
t~2 —> 73, in the second column. Similarly, in Tables IV and V show the obtained values of ^FA f° r ^ — 4 and 
J = 5, respectively. After carefully examining the previous tables, one may realize that there exists a fixed cost 
that depends on the number of broken constraints b that is b + 1 (at least, one call to W g WF is required before 
and after the broken constraint) and a variable cost that depends on the placement of the broken constraint. If the 
broken constraint is the last one the variable cost is 1. If it is the one before the last one, the variable cost is 2 
and so forth up to the case in which the broken constraint is the first energy causality constraint where the variable 
cost is J — 1. From this observation we can find E {C^^ FA } for a general J as 

,7-1 



b=0 
J-l 



J-l\ b + 1)+ fJ-2\(J-l)J 



b p ' \b-lj 2 



p^l-p)' 1 - 1 -" (41) 



= E V l 5 + + (42) 



b=0 



= ^ + l)(J-l)p+l=(y + ^-l)p+l, (43) 

where in (43), we have used that the mean of a binomial distribution with parameters n and p is np. Similarly, the 
variance of C^ S 1 ^ FA can be obtained through the variance of a binomial distribution as 

4- 1 ' 

2 



var{C£f4 FA } = ( ~ + 1 ) (J- l)p(l - p). (44) 



This concludes the analysis of the computational complexity of the algorithms. 



J=3 


Constraint 


^"HgWFA 


Probability 


✓ ✓ 

✓ X 
X ✓ 
X X 


1 

3 
4 
6 


(1-P) a 
(1 - p)p 
(1 - p)p 
P 2 


E {Cffi FA } = (1 - p) a + 7(1 - p)p + 6^ 



Table III 

Computational complexity of the FSA in the average case scenario (in terms of calls to H s WF) for J 



J=4 


Constraint 


WgWFA 


Probability 


✓ ✓ ✓ 


1 


(1-P) 3 


✓ ✓ X 

✓ X ✓ 
X ✓ ✓ 

✓ X X 
X ✓ X 
X X ✓ 
XXX 


3 
4 
5 
6 
7 
8 

10 


(l-pfp 
[l-pfp 
[l-pfp 
(1-P)P 2 
(1-P)P 2 
(l-p)p 2 


E {Cw^wfaI 


= {\-pY 


+ 12(1 - pfp + 21(1 - p)p l + 10p a 



Table IV 

Computational complexity of the FSA in the average case scenario (in terms of calls to H g WF) for J 



F. Properties of the reduction matrix 

The reduction matrix, Sk € xK , was introduced in [18] and is defined as: 



[S K ]i+(j-i)k,z = S ijz , {i, j, z} e [l,k] 
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J=5 


Constraint 


"Hjt WFA 


Probability 


✓ ✓ ✓ ✓ 

✓ ✓ ✓ X 

✓ ✓ X ✓ 

✓ X ✓ ✓ 
X ✓ ✓ ✓ 

✓ ✓XX 

✓ X ✓ X 
X ✓ ✓ X 

✓ X X ✓ 
X ✓ X ✓ 
X X ✓ ✓ 

✓ X X X 
X ✓ X X 
X X ✓ X 
X X X ✓ 
X X X X 


1 

3 

4 

5 

6 

6 

7 

8 

8 

9 

10 

10 

11 

12 

13 

15 


(i-rt 4 
(i -p) 3 p 

(i-p) 3 p 
(i -pfp 
{i -pfp 
(i -p?p 2 
{i-p)p 2 
{i -p) 2 p 2 
{i -pfp 2 
{i -p) 2 p 2 
{i -p) 2 p 2 
(i-p)p 3 
(i-p)p 3 
(i-p)p 3 
(i-p)p 3 
p 4 


E {C«*wfa} = (1 - Pf + 18(1 - pfp + 48(1 - p) V + 46(1 - + 15p 4 



Table V 

Computational complexity of the FSA in the average case scenario (in terms of calls to H g WF) for J = 5. 



Note that from the structure of S^, in each column there is only one entry different than zero and it is equal to 
one. For instance, the matrices for K = 2 and K = 3 are: 



(l \ 






The reduction matrix is designed so that 



and S? 



(l 







1 







v° 1 



S^vec(A) = diag(A) 



(46) 



for A £ $l KxK . In this appendix, we summarize some additional properties of the reduction matrix: 
Property 3. Multiplication properties: 
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• Let A G $l K xR , then the multiplication Sj~-A removes K 2 — K rows of A. 

• Let A G $i KxR , then the multiplication S K A adds K 2 — K rows of zeros to A. 

• Let A G $l RxK , then the multiplication AS]~- adds X 2 — A' columns of zeros to A. 
« Let A G SR^*^ 2 , then the multiplication ASk removes K 2 — K columns of A. 

Proof: The proof follows directly from the structure of the reduction matrix. ■ 
Property 4. Let A G U KxR , B G $t K xR , then S T K (A ® B) S K = A o B. 

Proo/: See [18, Lemma A.2]. ■ 
Property 5. S k Sk = lif • 

Proof: The proof directly follows from setting A = Ik and B = Ik in Property 4. ■ 
Property 6. Let A G $l KxK , then S^(A <g> I K )S K = Diag(diag(A)). 

Proof: The proof directly follows from setting B = Ik in Property 4. ■ 

Property 7. Let v G 5i A ', then Sj-(v ® I K ) = Diag(v). 

Proof: The Kronecker product expands the vector v in a K 2 x K matrix that stacks K diagonal matrices. 
Then, the multiplication by S K eliminates rows (see Property 3) so that the resulting matrix is Diag(v). ■ 

Property 8. Let A G '$t K * xK * be a diagonal matrix, then S^Sj-ASif = ASx 

Proof: From Property 3, S^-A removes rows from A. Then, the product by the left by Sk adds rows of zeros. 
As a result, S^S^A G %t R2xK ~ zeroes K 2 — K rows of A. Finally, the product with from the right removes 
K 2 — K columns. As A is diagonal, the entries that are modified by multiplying from the left by S^Sj- are later 
removed by multiplying from the right by Sk- Therefore, Sa'S^-ASk is equal than AS^, which directly removes 
the columns. ■ 
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